THERMOHALINE MIXING: DOES IT REALLY GOVERN 
THE ATMOSPHERIC CHEMICAL COMPOSITION OF 
LOW-MASS RED GIANTS? 



Pavel A. Denissenkov 1 , and William J. Merryfield 2 
pavel . denisenkov@gmail . com 
bill .merryf ield@ec .gc . ca 

ABSTRACT 

First results of our 3D numerical simulations of thermohaline convection 
driven by 3 He burning in a low-mass RGB star at the bump luminosity are 
presented. They confirm our previous conclusion that this convection has a mix- 
ing rate which is a factor of 50 lower than the observationally constrained rate 
of RGB extra-mixing. It is also shown that the large-scale instabilities of salt- 
fingering mean field (those of the Boussinesq and advection-diffusion equations 
averaged over length and time scales of many salt fingers), which have been ob- 
served to increase the rate of oceanic thermohaline mixing up to one order of 
magnitude, do not enhance the RGB thermohaline mixing. We speculate on pos- 
sible alternative solutions of the problem of RGB extra-mixing, among which the 
most promising one that is related to thermohaline mixing is going to take advan- 
tage of the shifting of salt-finger spectrum towards larger diameters by toroidal 
magnetic field. 

Subject headings: stars: abundances — stars: evolution — stars: interiors 

1. Introduction 

Red giants with M < 1.5 M & are known to experience extra- mixing in their convectively 
stable radiative zones that separate the hydrogen burning shell (HBS) from the bottom 
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of the convective envelope (BCE). It operates during their first ascent along the red giant 
branch (RGB) above the bump luminosity, after the HBS has crossed a chemical composition 
discontinuity left behind by the BCE at the end of the first dredge-up. This causes the 
atmospheric 12 C/ 13 C ratio and carbon abundance to resume their declines with an increase 
in luminosity. 

The most promising physical m echanism of RGB extra-mixing is thermohaline con- 
vection (ICharbonnel fc Zahnll2007al ). It is driven by a double-diffusive instability that oc- 
curs when a scalar component that stabilizes a density stratification, e.g. the temperature 
T, diffuses away faster than a destabilizin g component, e.g. the mean molecular weight 
[i. lEggleton. Dearborn, fc Lattanziol (120061 ) brought our attention to the fact that the re- 
action 3 He( 3 He, 2p) 4 He produces a local depression of /i, A/i ~ — 10~ 4 , in the HBS tail 
above the bump luminosity. As a result, fi is increasing with the radius, which turns the \x 
profile into the destabilizing component. Moreover, the radiative thermal diffusivity K ex- 
ceeds the molecular diffusivity z/ mo i by many orders of magnitude (the inverse Lewis number 
t = u mo \/K ~ 10~ 6 in the region of \i depression), which should lead to the double-diffusive 
instability. In the ocean, a similar instability develops where warm salty water overlies cold 
fresh water. The oceanic thermohaline mixing usually takes the form of vertically elongated 
fluid p arcels of ris ing fresh and sinking salty water adjacent to each other, called "salt fin- 
gers" JStern ]l960|) . It has extensively been studied both experimentally and theoretically 
(e.g., 



Kunze 



2003 



The efficiency of thermohaline mixing in stellar cases can b e estimated o nly theoreti- 
cally, e.g. via a linear stability analysis of its governing equations (IUlrichlll972l ) or via direct 
numerical simulations. Denissenkov (2010, hereafter Paper I) has recently conducted a com- 
parative study of thermohaline convection in the oceanic and RGB cases using both a linear 
stability analysis and 2D numerical simulations. Unfortunately, the linear theory yields a 
thermohaline diffusion coefficient that is proportional to the square of the (unknown) aspect 
ratio, a = l/d, of a finger, where I and d are the finger's length and diameter. Applica- 
tions of this diffusion coefficient to model the evolutionary decline of carbon abundance in 
low-mas s red giants demand t hat a > 7 (Paper I), which is in agreement with the original 
result of ICharbonnel &: Zahnl (j2007al ). However, our 2D numerical simulations of 3 He-driven 
thermohaline convection in the vicinity of a fi depression yield a mixing rate equivalent to 
that approximated by the linear-theory diffusion coefficient with a < 1, i.e. it is a factor of 
~ l/50 th of its observationally constrained value. The difference in the finger aspect ratios 
between the oceanic (a > 1) and RGB (a < 1) cases is most likely determined by the very 
low viscosity u, or the Prandtl number Pr = v/K~ 10 -6 , in the RGB case. This facilitates 
the developm ent of secondary shear instabilities that strongly limit the growth of salt fingers 
( IRadkol 120101 . and Paper I). On the other hand, the diffusion coefficient for thermohaline 
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convection with a > 7 has been shown to reproduce many of the RGB extra-mixing observa- 
tional patterns at all metallicities so well that we had to postpone further discussion of these 
findings until we were able to confirm or refute the results of our 2D numerical simulations 
using 3D computations. 

In this Letter, we present first results of our 3D numerical simulations of thermohaline 
convection for the oceanic and RGB cases. It turns out that they do not differ much from 
the 2D results reported in Paper I. We also consider a possibility, mentioned in Paper I, that 
RGB thermohaline mixing may be enhanced by large-scale mean-field instabilities, such as 
the collective, 7-, and intrusive ones, that are known to occur in the ocean. Finally, we 
speculate on possible alternative solutions of the problem of RGB extra-mixing. 



2. 3D Numerical Simulations of Thermohaline Convection 

The basic equations and numerical techniques that are employed in our 3D simula- 
tions of thermohaline convection for th e oceanic and RGB cases are described in detail by 



Gargett. Merryfield. fa Hollowayl (120031 ). They are very similar to those used in Paper I. The 
main goal of Gargett et al. was to study differential scalar diffusion in 3D stratified turbu- 
lence. Therefore, they considered a density stratification stable against the double-diffusive 
instability and used a specified turbulent velocity field as an initial condition. To apply their 
method and computer code for our purposes, we have made the slower diffusing component, 
which represents salinity S in the ocean and /1 in the current investigation, destabilizing by 
changing the sign of the vertical velocity term in their equation (16). We have also reduced 
the amplitude of the initial velocity field, so that it serves as a mild perturbation to initi- 
ate the double-diffusive instability. Our available computational resources have allowed us 
to run the code within reasonable time intervals with a spatial resolution of 320 3 . Simple 
estimates of the relevant turbulent kinetic energy dissipation rate and the Kolmogorov and 
Batchelor length scales show that this resolution is close to marginal for our runs, however, 
not fully resolving the sl ower- diffusing comp onent has been shown not to have a large effect 



on the computed fluxes (ITraxler et al.ll2010l ). The 3D simulations have been carried out for 
the same number of fastest growing fingers in the computational domain and parameter sets 
from Table 1 in Paper I that were used in our 2D computations. 

The blue and red solid curves in the upper left panel of Fig. [TJshow transitions of the ratio 
Ds/kr to its equilibrium values in our 2D and 3D simulations of the oceanic thermohaline 
convection. Here, Ds and kx are the turbulent salt and microscopic thermal diffusivities. 
The 2D simulations were performed with the resolution 1024 2 (Paper I). Our estimated 3D 
equilibrium value of Ds/k? exceeds its corresponding 2D value merely by a factor of 1.5 
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(compare the blue and red dashed lines). It is i n good agreement with the result obtained 
using a better resolution by iTraxler et al.l (120 lOf ) that we have interpolated for the value of 
density ratio R p = aAT/f3AS = 1.6 (a and (3 are the coefficients of thermal expansion and 
haline contraction) from their Table 1 (star symbol). The upper right panel makes a similar 
comparison for the RGB thermohaline convection. In this case, we do not find any difference 
between the 2D and 3D equilibrium values of the ratio D^/K (the blue and red dashed lines 
coincide). Therefore, we conclude that our results reported in Paper I seem to be correct. In 
particular, RGB thermohaline convection driven by 3 He burning yields a turbulent mixing 
rate that is by a factor of 50 lower than the value required by observations. 



3. Mean-Field Instabilities 



Observations, laboratory experiments, and direct numerical simulations show that small- 
scale salt-fingering convection in the ocean can lead to instabilities and the formation of dy- 
namical structures on much larger scales. These include the so-called "collective instability" 
leading to the generation of internal gravity waves, the intrusive instability developing when 
fluid is stratified both vertically and horizontally (in stars, such a situation may result from 
the rotational distortion of level surfaces), and the s o-called "7-instability" that produces 



thermohaline staircases (e.g., see ITraxler et al.l I201CH . and references therein). The latter 



structures consist of thick well-mixed convective layers separated by thin salt-fingering inter- 
faces. The staircases are of special interest to us because the y have been observed to enhance 
vertical mixing in the ocean by up to an order of magnitude (ISchmitt et al.ll2005l ). Moreover, 
our preliminary estimates in Paper I indicate that the turbulent flux ratio 7 = Fx/F p is a 
decreasing fu nction o f R p i n the RGB case, which is a sufficient condition for the 7-instability, 
according to iRadkol (120031 ) . 



The large-scale instabilities are studied by averaging the Boussinesq and advection- 
diffusion equations for T and S (//) over length and time scales of many salt fingers. As 
a result, a system of mean- field equations is obtained in which the small-s cale turbulent 



veloc ity is incorporated into the mean heat and salinity fluxes Ft and Fg ( ITraxler et al. 



20101 ) . The fluxes are related to one another and to the large-scale temperature gradient 
via the flux ratio 7 = Ft/Fs and the Nusselt number Nu, which are assumed to depend 
only on the local value of the density ratio R p ; specification of these dependences closes the 
mean-field equations. Large-scale changes in R p thus result in modulations of F T and F s via 
the dependences of 7 and Nu on R p that can either amplify or damp initial perturbations. 

A linear stability analysis of the mean-field equations yields a third-order dispersion 
relationship for the growth rates of all the three large-scale instabilities mentioned above 
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(equation 2.13 of iTraxler et al.l |2010| ). Coefficients of the cubic equation depend on the 
Prandtl and inverse Lewis numbers as well as on Nu, 7 and their derivatives with respect 
to R p . The coefficients are also functions of the horizontal and vertical wavenumbers, I 
and k. To ca l culate them, we first need to estimate the dependences of Nu and 7 on R p . 



Traxler et al.l (120101 ) estimated those dependences via a body of 3D numerical simulations of 



the oceanic salt-fingering convection. We have used data from their Table 2 to reproduce, for 
a comparison, their "flower plot" depicting the growth rates of fastest growing instabilities 
at R p = 1.5 in the absence of lateral gradients (the lower left panel of Fig. [T]). The bulb 
of the flower corresponds to the salt-fingering mode that turns out to be a solution of the 
general dispersion equation as well. The leaves of the flower represent the gravity waves 
of the collective instability, whereas the 7-instability modes are the ones dominating in 
the region beneath the leaves. To obtain estimates of 7(-R p ) and Nu(R p ) for the RGB 
case (Fig. [2]), we have performed additional 2D numerical simulations of the 3 He-driven 
thermohaline convection for the same parameter set that was used in Paper I but for a 
number of different values of R p . The resulting flower plot is shown in the lower right panel 
of Fig. [H We see that in the RGB case only the salt-fingering mode is unstable in spite of 
the flux ratio 7 being a decreasing function of R p , as in the oceanic case. This difference 
is explained by the fact that the Nusselt number is extremely small, Nu ~ 10~ 6 , in the 
RGB case. Because of this, it is necessary to include horizon tal diffu s ive flu xes in the mean- 
fiel d stability ana l ysis — something which was not done by iRadkol (120031 ) but was treated 
by ITraxler et al.l (j2010l ). In the RGB case, the inclusion of reasonable lateral gradients 
proportional to their corresponding vertical gradients multiplied by the ratio of centrifugal 
to gravitational acceleration only produces an asymmetry of the bulb with respect to the 
vertical axis. 



4. Possible Alternative Solutions 

Our 2D and 3D numerical simulations of thermohaline convection driven by 3 He burning 
have demonstrated that it is unlikely to be the sole mechanism of RGB extra-mixing. Its 
corresponding mixing rate in the vicinity of a u depression, si2x 10~ 3 K (the upper right 
panel of Fig. [I]), is too low compared to the observationally constrained value of (-D A1 ) bs ~ 
0.1 K (Paper I). Therefore, for its efficiency to be consistent with observations, the RGB 
thermohaline mixing has to be either assisted or modified by other processes. In Paper I, 
we have shown that a turbulent viscosity u t exceeding the microscopic one by a factor of 
10 4 could enhance the rate of 3 He-driven thermohaline mixing up to the observed value, 
provided that it was not itself associated with strong turbulent mixing. The enhancement 
of thermohaline mixing in this case is caused by the fact that the higher viscosity suppresses 
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the development of secondary shear instabilities that limit the growth of salt fingers. On 
the other hand, if the increase of viscosity is accompanied by a similar increase of the rate 
of turbulent mixing, the latter will reduce salt-finger buoyancy and vertical transports by 
smoothing out the chemical composition contrast between rising and sinking fluid parcels. 

The additional source of turbulence in the radiative zones of RGB stars could arise from 
the differential rotation that is produced by the mass inflow that feeds the HBS (Paper I). We 
have employed the COMSOL Multiphysics software package to solve the angular momentum 
transport equation for the radiative zone of the low-metallicity bump-luminosity RGB model 
from Paper I in the presence of mass inflow and turbulent shear mixing. Their respective 
rate r and diffusion coefficient v w are given by equ ation (23) from Paper I and equatio n 



(3) with the uncertainty factor f v from the paper of iDenissenkov. Chaboyer. fc Li I ( 120061 ). 



The angular velocity at the BCE was assumed to h ave the same value, ^rcf, = 6.3 x 10 
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rads -1 , as in the M2 bump- luminosity RGB model of lPalacios et al.l (120061 ) . Simple estimates 
show that for / v > 1 the transport of angular momentum by the rotation-induced turbulent 
diffusion dominates over its advection by meridional circulation, therefore, and for the sake 
of simplicity, the latter was neglected. Our computations show that the competition between 
the terms proportional to r and v N leads to a stationary solution with u v ~ 3 x 10~ 4 K in the 
vicinity of a fi depression that very weakly depends on both f v and ^bce, although the final 
rotational shear is a decreasing function of f v . This value exceeds the microscopic viscosity 
by a factor of 10 2 , which is not enough to enhance to its observed value (Paper I). It 
is also difficult to understand why the rotation-induced small-scale turbulence should not 
produce a mixing rate of the same magnitude as u v , in which case will not be enhanced 
at all (Paper I). Therefore, we do not think that the proposed viscosity enhancement is a 
possible solution of the problem. 

Our angular momentum transport computations confirm the conclusion that radiative 
zones of low-mass RGB stars above the bump luminosity should possess strong differen- 
tial rotation. This may lead to the generation of high-amplitude toroidal magnetic fields, 
especially close to the HBS, provided that a weak poloidal field is present in the radia- 
tive zone. Under these circumsta nces, a possible alternative solution of the proble m could 



be magneto-thermohaline mixing ( iDenissenkov. Pinsonneault. &: MacGregor 1120091 ) . We in- 



pi 

2> 



tend to investigate this possibility using 2D numerical simulations in one of our forthcoming 
papers. 

Toroidal magnetic fields in radiative zones may also filter out salt fingers with small 
diameters, so that the maximum growth rate is shifted towards thicker fingers. Since the 
finger growth rate decreases with an increase in diameter as o oc dr 2 (e.g., Paper I), the 
square of the vertical velocity shear between the axes of neighbouring rising and sinking 
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fingers, ~ (al/d) , is proportional to a d . This par ameter is likely to control the de- 



velopment of secondary shear instabilities (IKunzd 119871 ); consequentl y, thicker fingers will 



probab ly reach higher aspect ratios before they are destroyed. Following ICharbonnel fc Zahn 



(l2007bl ). we have complemented the dispersion relationship (13) from Paper I (note that the 
term uhrksk 6 was erroneously omitted) with terms that take into account the Lorentz force 
associated with a specified toroidal magnetic field B . Its solutions are shown in Fig. [3] for 
the field strengths B = 0, 10, and 100 G (the solid blue curves from left to right). In the 
case of Bq = 100 G, the fastest growing fingers transverse to the field have a diameter that is 
nearly 40 times as large as that of fingers in the absence of a magnetic field. This diameter 
is also almost four times larger than the diameter of fastest growing fingers in the case of a 
viscosity that is amplified by the factor of 10 4 (the red solid curve). However the influence 
of such a field on the planform of the fastest growing fingers still needs to be studied. In our 
next study, we will perform 2D numerical simulations of the RGB thermohaline convection 
in the presence of a toroidal magnetic field to find out if the magnetic shift of the salt-finger 
growth-rate spectrum towards larger diameters can really enhance its associated mixing rate 
up to the empirically constrained value. 
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Fig. 1. — Upper panels compare the ratios of thermohaline diffusion coefficient to thermal 
diffusivity obtained in our 3D (red curves) and 2D (blue curves) numerical simulations for 
the oceanic (left) and RGB (right) cases. Black star symbol shows a value interpolate d 
for the density ratio R p = 1.6 using the higher resolution 3D data of iTraxler et al.l ( 120101 ). 
Lower panels depict the decimal logarithms of growth rates (s -1 ) of the fastest growing large- 
scale instabilities of salt-fingering mean-field as functions of the horizontal I and vertical k 
wavenumbers, both measured in cm -1 , for the oceanic (left) and RGB (right) cases. 




Fig. 2. — The inverse flux ratio, 7 1 = F^/Ft, and Nusselt number as functions of density 
ratio estimated from our 2D numerical simulations of the RGB thermohaline convection. 
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Fig. 3. — The RGB salt-finger growth-rate spectra computed for the same parameter set 
as in Table 1 from Paper I under the assumptions that the toroidal magnetic field has the 
strengths B = 0, 10, and 100 G (the blue curves from left to right), and for B = G but 
with the viscosity amplified by the factor of 10 4 (the red curve). 



